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Abstract 

We study the efficiency of synchronization in ensembles of identical coupled stochastic oscillator systems. 
By deriving a chemical Langevin equation, we measure the rate at which the systems synchronize. The 
rate at which the difference in the Hilbert phases of the systemsevolve provides a suitable order parameter, 
and a 2-dimensional recurrence plot further facilitates the analysis of stochastic synchrony. We find that 
a global mean-field coupling effects the most rapid approach to global synchrony, and that when the 
number of “information carrying” molecular species increases, the rate of synchrony increases. The 
Langevin analysis is complemented by numerical simulations. 

Keywords: Cell signaling, synchronization, diffusive coupling, mean-field coupling, Neurospora 
crassa. 

Introduction 

The nature of synchrony in stochastic dynamical systems has been a subject of interest given the ubiquity 
of both stochasticity and synchrony (or more generally, strong temporal correlativity) in a variety of 
natural systems m- The dynamics of biological systems, particularly at the cellular or subcellular 
level, are known to be subject to large fluctuations, and it is clear that diverse processes need to be 
fine-tuned temporally in order that any biological “system” is able to function. Other examples can be 
found in areas ranging from neuroscience to the nature of financial markets, and indeed, in other instances 
when a systems approach is applicable [IHS]. 

The coupling of autonomous dynamical systems can be effected through a variety of different mecha¬ 
nisms such as subjecting them to a common driving signal [7], or through diffusion [HHin] or through a 
mean-field [am]. One of the most widely studied effects of such coupling is the emergence of new collec¬ 
tive behaviour, for instance synchronization m- As has been established through extensive work in the 
past decades, synchronization comes in a number of flavours: phase [12) . diffusive [13) . frequency [14) . 
delay [15) , complete [7] , generalized [16) synchronization, based on which variables are being synchronized 
and how the synchrony is manifest. Different measures or order parameters to gauge the rate of syn¬ 
chrony under the various scenarios have been devised: for instance the phase difference A(f> as a function 
of time [12], Lyapunov exponents permutation entropy phase locking value [18] etc. Several of 
these scenarios apply in the presence of noise, both external and internal [TTHT^ISO] . 

In the past few years, the study of stochastic synchronization, namely the emergence of temporal 
correlations between stochastic dynamical systems has been a major area of study [3]. The approach to 
phase synchrony in such cases depends both on the nature of the coupling as well on details relating to 
time delay or coupling topology, and this issue forms the central question of the present paper: How does 
synchronization efficiency depend on these different factors? 
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This is a central issue in the related context of how the dynamics in a group of systems become 
correlated. To give a biological instance, circadian rhythms in bacterial organisms such as Neurospora 
crassa are globally correlated through environmental fluctuations. A number of biological rhythms are 
based on genetic processes that originate in negative feedback circuits to generate a periodic 

expression of specific genes within a cell [2T1[^. Clock proteins like FRQ in N. crassa, or PER and 
TIM in D. melanogaster are responsible for the genetic regulation of clock genes [261127] . Environmental 
fluctuations are known to be a means of correlating a group of cells or unicellular bacteria, most notably 
via the mechanism of quorum sensing [28] . Synchrony appears to be essential for information processing [1] 
and is the consequence of celEto-cell communication via specific coupling mechanisms [29] . 

In order to focus on the measurement of the rate of synchrony for various coupling mechanisms 
and to identify the factors on which this rate depends, we study coupled stochastic oscillators within 
the “chemical Langevin” formalism. This is presented in Section II for different coupling mechanisms. 
In Section III we examine a circadian oscillator model which we study in detail both numerically and 
analytically. Our results are summarized in Section IV. 

Chemical Langevin eqnation formulism of coupled oscillators 

The random interaction of molecules in a well stirred mesoscopic system leads the dynamics of the 
variables in the system to noise-driven stochastic process [301134] . Consider the state of the system at 
any instant of time t is defined by a configurational state vector, C{t) = [Xi{t),X 2 {t),...,XM{t)Y', where 
N distinct molecular species are interacting via M elementary reaction channels in the system of the 
following type, 


oti^Xi H- aNf^Xpf -4- (5\^Xi H- Pn^iX]^ (1) 

where, /r is the reaction number index: fj, = 1,N, cxi^ and are co-efficients to define stoichiometric 
matrix and is the /xth macroscopic rate of reaction. The system evolves with various 

random reaction fired at random interval of time with decay or/and creation of molecular species at any 
reaction event [301131] which leads to the change in configurational state of the system. This allows to 
define a configurational probability P(C;t) as the probability to get this state change in an interval of 
time [t,t + At]. Then the time evolution of configurational probability P{C;t) obeys Chemical Master 
equation (CME) [501135] . The CME in general provides detail mesoscopic description of chemical kinetics, 
but it is very difficult to solve for complex systems [5T] . 

The chemical Langevin equation (CLE) formalism is one method to approximate CME to simpler 
continuous Markov type equations by keeping conditions which are applicable in natural systems [32] . 
and the accuracy of this CLE is found to be more than those of other formalisms such as linear noise 
approximation [36] . The approximation can be done by allowing to define a function A{C,At) as the 
number of a particular reaction fired during an interval of time [t, t At] with At)0. This is followed by 
excellent approximations by imposing two conditions, firstly, imposing small At limit such that the values 
of propensity functions w(C) of the reactions remain constant during [t,t -|- At], and secondly imposing 
large At limit which in turn leads to a;(C(t))At))l. These two conditions allow A to approximate to 
statistically independent Poisson random variable and then the Poisson random variable is replaced by 
normal variable with the same mean and variance. Both the conditions are true in natural practice 
for large population limit. Then linearizing the normal variable, and defining macroscopic molecular 
concentration vector {x} = where V is the systems size, we have general CLE, 

^ = E[u;(x(t)), H + G[a;(x(t)), z., E] (2) 

where, F = VijUJi[x{t)\ is the macroscopic contribution term and G = 
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the stochastic contribution term to the dynamics, = limdt^oNiiQ, l)l\/dt is uncorrelated, statistically 
independent random noise parameters which satisfy = SijS(t — t'). 

Now consider two identical interacting stochastic systems with configurations U{t)=[xi{t)^X2{t)^...^XM(ty\^ 
and U' {t)=[x'i{t),x' 2 {t),...,x'j^{t)Y', where their dynamics are given by CLE of the type equation (I2|), 
x= U{x) and x = U'{x') respectively. The interaction of the two stochastic systems can be allowed by 
choosing one or more ’’coupler” molecular species and introducing one of the various coupling mech¬ 
anisms, namely, direct (master-slave) [7], diffusive [S], mean-field m coupling etc via y. We then 
construct a larger 2N-dimensional stochastic system, H{z) whose dynamics is given by, z = H{z) where, 
z = (xi,..., xn, x(l, ..., x(^). Then the system is devided into sub-systems, U{x) and U'{x') consisting of 
independent reaction sets with corresponding dynamics and a third arbitrary reducible sub-system S{y) 
in the same configurational space H{z) formed by m extra reaction channels introduced via 2n coupling 
molecular species y = (cci+i,..., cci+n, • ■ •; ■ If there is no coupling between the two systems, 

S' —>■ 0, otherwise S is finite and it depends on various factors and parameters such as directionality of 
coupling, coupling constants etc. This S associates with signal or signals common or diffused between 
the two sub-systems derived from the logical operations or extra reaction channels responsible for the 
coupling. This leads to the following dynamics of the S', U and U' sub-systems given by. 


dy{t) 

dt 

dx{t) 

dt 

dx'{t) 

dt 


^(y,i^,e,V) + S,(y,V,e,0 


f7(x,i/,^, V), (xT^y) 

f7'(x',z/,^, V) (x^y) 


(3) 

(4) 

(5) 


where, U = [H{xi),..., H{xn)]^ and U' = [H{x'l),..., H{x'j^)]'^ are of the form of equation ([ 2 |) for 
uncoupled casees. The last term Sy(y, V, e,f) is the coupling term added to the corresponding coupling 
variables, e is the coupling constants between the two subsystems. The functional form of Sy depends 
on how the two systems coupled and type of coupling mechanism, for example: (i) Direct coupling : if 
Xj = x'j, (j = i + 1,... ,i + n)^ then Sy ^ 0 and the other variables will achieve synchronization [7], (ii) 
Diffusive coupling : if the diffusive reactions are incorporated to couple the two systems, then Sy is 
finite with those number of diffusive reactions (uni/bi directional) to form Sy and synchronization can be 
achieved among the rest of the variables for sufficiently large values of e, (iii) Mean — field coupling : if 
mean information of the two systems is allowed to use for signal transduction between the two systems by 
constructing an arbitrary molecular species, rjj = = 2) and allowed to diffuse to interact 

with corresponding molecular species, xj and x' in the two sub-systems via diffusive reactions, then 
synchronization will be exhibited among the other variables, and so on. 

We can generalize this coupling method for L identical stochastic systems by constructing a large 
stochastic system H = ... ,U^[N,M),Sy{Mii)], where is the total number 

of extra reaction channels allowed depending upon the type of coupling and the way how they couple 
among them. If we look for steady state solutions of CLE (3)-(5) by applying conditions, x= 0, i = 0 
and y= 0, then we obtain Sy{M}i) = 0. 


Synchronization efficiency: diffusive coupling 

We first consider diffusive coupling mechanism between two stochastic systems to understand the rate syn¬ 
chronization among the coupled systems. If y = (xq, x'^) is taken to be couplers which can diffuse in and 
out of the systems (bidirectional), this coupling between any two stochastic systems U = [xi, 0 : 2 ,..., 
and U' = [xi,X 2 t ...,x'j^\~^ can be achieved by incorporating two extra reaction channels; Xa A x'^] 

x'^^ Xa. When the coupling constants c and d are strong enough, the dynamics of other variables of 
the systems i.e. {{xi,x!i^-,i = 1,2,..., N,i ^ a} will achieve synchronization. In this case H will have 2N 
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variables, 2M + 2 reaction channels: M reaction channels for each sub-systems U and U', and 2 forS'y(2). 
We take c = c' for simple case and the functional forms of two Sy = , Sx '^) for two variables are given 

by the following CLEs derived from equation ([2]) for the two reactions, 


Sxa.iXa,x'^,V,C,^) 
= c[x'^ - Xa\ -b 


Sx'^{x'a,Xa,V,C,Cj 



= c[Xa 



( 6 ) 


Substituting these functions Si and S 2 in equations (3)-(5) we have CLE for the two coupled stochastic 
systems. Now we look solution of the CLEs for stability condition. Since x= U{x), x = U'{x) and H{y) 
is a part of U and U' for CLE (3)-(5) and (6), we have the stability conditions, x= 0, x= 0. This 
conditions give us Sx^ = 0 and Sx'^ = 0 respectively and after solving these two equations, we get the 
following result. 


1 


{Xa, , Xq^') 


(7) 


where, Ks{xa,x'^)= 



For the sake of simplicity we take 




VB or Cs 




where i? is a constant and then taking Xa x'^, we have, K^ixa^ x'^) ^ So for bidirectional diffusive 
coupling of single molecular species, the coupling constant is related to system’s size by, c ~ 

Now we consider two molecular species Xa and Xjs as ’’couplers” and are allowed to diffuse among the 
two systems U and IB in the similar fashion discussed above by constructing Sy{y) = {xa:Xp,x'^,x'pY' 

c c c" c" 

by introducing four extra reaction channels: x^ —>■ x'^\ x'^ x^] xp —>■ x'p —>■ xp giving S = 

{Sx^, Sx'^, Sxp, Sx'^)- We then substitute these forms of Sy in the equation (3)-(5), and the equations 
become CLE for the two coupled stochastic systems. We then look for steady state solutions by applying 
steady state conditions which give all four 5" = 0. Solving the equations for Xa and x'^ and keeping c = c' 
we have the relation between c and V for is given by equation ©• Similarly keeping c" = c'", the 
relation between c" and V for xp can be obtained by solving the equations 5 = 0 for xp and x'p which 
leas to the following equation, 


-^=K{xp,x'p) (8) 

where, we take ~ ~ (Cd ^ ^'e ^ V^) and then taking xp ^ x'p, we have, A',,{xp,x'p) ^ 

such that c" ^ y. 

The rate constant of a chemical reaction is the product of number of collisions per unit time and the 
probability that any given collision in the reaction takes place [38] . If the two sub-systems are coupled via 
diffusion of two molecular species with two different coupling constants c and c", the equivalent coupling 
constant c' of the two diffusing rate constants will follow ’’parallel resistance law” hypothesis [39l - l^ 
given by ^ i -b ^ 7 . The hypothesis can be generalized for L different rate constants given by vector 
c = (ci, C 2 ,..., cl)^ to obtain equivalent rate constant - = —-b — -b----b—. From the equations dTl) 

C Cx C2 

and ([8|), the equivalent coupling constant c is solved using this hypothesis, 


1 

Xa+Xfi 

where, we take B = B' for the sake of simplicity and since they are random numbers. 



( 9 ) 
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Synchronization efficiency: mean-field coupling 


We now consider E identical subsystems, and define an arbitrary hypothetical mean field molecular 
species, ria{t) = which is the average information carried by K subsystems and allowed 

to signal transduction via molecular species Xa in the ensemble. Depending on the number A and the 
topology of the ensemble of E sub-systems, the mean-field coupling could be nearest neighbour {K =2 
(1-dimension), 4(2-dimension), 6(3-dimension)), next to nearest neighbours and global [K = E). Any 
two oscillators Uj and C/,+i in the ensemble of E sub-system can be coupled via mean-field by allowing 

and to diffuse between the two subsystems via two extra reactions, 77^+1 

Thus we can construct 5'(a;^,because 77 ^ = rjl^ix^^) and 77 ^+^ = This mechanism is 

similar to that of diffusive kind but exchange of mean information takes place. Proceeding in the same 
way as we did in diffusive coupling case, we can get the set of CLE by substituting the calculated S in 
(3)-(5). Taking e = e' the relation between e and V is obtained as. 


1 1 D 

An{ri,r]')V E 477 


( 10 ) 


where, A„{r], ri')= 


and (e^ ~ ^ VD). 


We then couple the two oscillators and C/ 2 I -1 ensemble with mean-field coupling mechanism 

via two molecular species and x^^^ by defining two corresponding hypothetical mean-field molecular 

species ria{t) = and 7 /i^(t) = respectively. These rja and ipp are allowed 

to diffuse in the sub-systems and interact with the local molecular species Xa and xp in the respective 


sub-systems. This can be done by defining four extra reactions: 77 ^ 


’ 'la 




j+^. 


^ This leads us to construct 5'(x^,x^,a:^"'’^) because rja and ipp depend on Xa and xp 
respectively. Then substituting S forms in equations (3)-(5) we get CLE of the coupling mechanism. 
Then putting e = e' for ija and e" = e'" for 'ipp, we solve the CLE for steady state conditions and using 
the” parallel resistance law” hypothesis for e and e", we get the equivalent rate constant e' as in the 
following, 


e 


/ 


D 1 
4V 7] + ip 


( 11 ) 


As we see from the above equation that the functional form of the rate constants in both mean field and 
diffusive couplings are similar. 


Comparison of the coupling mechanisms 

The comparison of the coupling constants in all the four coupling mechanisms is done so that one can 
compare the rate of synchrony achieved between the two subsystems. From equations © and we 
found ^ ~ ^'^+^17 ^ 2 ^ giving rise c)c'. Again from equation (0 and (fTUj) . we could arrive at ^ ~ yT^(l 

with B ^ D such that c'(e. Similarly, from equation (ITUl) . (fTT|) . ([5]) and (fTUl) . we can show that jr ~ ^^)1 
and ^ ~ respectively giving rise e)e' and c' > e'. Thus summarizing the possible relations, we 

have, 

c > e) c' > e' ( 12 ) 

This indicates that the synchronization of the oscillators occurs at largest value of coupling constant in 
the case of single molecule diffusive coupling and at the smallest value of it in the case of double molecule 
mean-field coupling. In another words, the rate of information processing from one oscillator to another 
is fastest in the case of double molecule mean-field coupling mechanism and the rate of this information 
processing in ascending order is given by equation (Ha. 
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Figure 1. Biochemical reaction network of circadian rhythm in N. crassa due to Gonze and 
Goldbeter [44] . 
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Figure 2. Coupling mechanisms we considered in our simulation: (A) Nearest neighbour single 
molecular species (M is coupling molecular species) bidirectional diffusive coupling, (B) Nearest 
neighbour single molecular species {ti = ^ Sti coupling molecular species) bidirectional 

mean-field coupling, (C) Nearest neighbour double molecular species (M and Fc are coupling molecular 
species) bidirectional diffusive coupling and (D) Nearest neighbour double molecular species 
iv = i Eh and are coupling molecular species) bidirectional mean-field 

coupling. 
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Figure 3. Plots of nuclear protein Fjv(t) for 10 oscillators out of 50 oscillators as a function of time, 
for V = 100 and e = 1.4 showing synchronized and desynchronized regimes. Phase differences of pairs of 
oscillators, i,j = {(1,2), (1,3),(1,4), (1,5), (1,6), (1,7), (1,8), (1,9), (1,10), (1,50)} as a function of 
time i.e. phase plot, phase lock values of the corresponding set pairs of oscillators and recurrence plot 
are also shown in second, third and fourth panels respectively. 


Measurement of rate of synchrony 

It has been pointed out that the identification of phase synchrony of any two identical systems can be 
done qualitatively by the measurement of the time evolution of instantaneous phase difference of the two 
systems [2J[9l[20l|42] . It is possible if one defines an instantaneous phase of an arbitrary signal x{t) via 
Hilbert transform [20] 


x{t) = -P.V. [ (13) 

J -oo t ~ T 

where P.V. denotes the Cauchy principal value. Then, one can determine an instantaneous ’’phase” 
and an instantaneous ’’amplitude” A{t) of the given signal through the relation, = x(t) + ix{t). 

Given two signals of two systems Xi(t) and Xj{t), one can therefore obtain instantaneous phases 
and <(>j(t); phase synchronization is then the condition that Acjjij = mc/ii — n(j)j is constant, where m and 
n are integers. Of most interest are the cases A(j)ij = 0 or tt, namely the cases of in-phase or anti-phase, 
but other temporal arrangements may also occur. 

The phase synchronization of the two identical systems can also be identified by doing synchronization 
manifold recurrence plot of the variable of one system, say x with the corresponding variable x' of the other 
system simultaneously on two dimensional cartesian co-ordinate system [Tj. The rate of synchronization 
between the two systems can be determined by the rate of concentration of the points in the plot along 
the diagonal. If the two systems are uncoupled, then the points on the plot will scatter away randomly 
from the diagonal. 

The rate of phase synchrony of the two systems can be estimated quantitatively by measuring the 
’’phase locking value” (PLV) of the two signals of the two systems [18]. The phase locking value, which 
is used to quantify the degree of synchrony, characterizes the stability of phase differences between the 
phases 4>i (t) and 4>j (t) of two signals Xi (t) and Xj (t) of ith and jth systems and can be defined within a 
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Figure 4. Plots showing transition from desynchronized to synchronized regime for different coupling 
mechanisms at F = 100 and c(c or d or e, or e') = 0.6: (1) Diffusive coupling (single molecule): M is 
taken as coupling molecule. Phase plots, phase locking values {PLV) and recurrence plots are shown in 
panels (a), (b) and (c). (2) Diffusive coupling (double molecule): M and Fc are taken to coupling 
molecules. The respective plots are as in (1) are shown in panels (d), (e) and (f). (3) Mean-Field 
coupling (single molecule): M is taken as coupling molecule and mean field for 30 oscillators is taken. 
The respective plots as in (1) are shown in panels (g), (h) and (i). (4) Mean-Field coupling (double 
molecule): M and Fq are taken as coupling molecules and the corresponding plots are shown in panels 
(j), (k) and (1). 
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Figure 5. Similar plots as in Fig. 4 are shown but as a function of coupling parameters; Nearest 
neighbour (1) Diffusive coupling (single molecule) in (a), (b) and (c); (2) Diffusive coupling (double 
molecule) in (d), (e) and (f); (3) Mean-Field coupling (single molecule) in (g), (h) and (i); (4) 
Mean-Field coupling (double molecule) in (j), (k) and (1). 
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Figure 6. Phase plots {Acj) vs t), PLV and recurrence plots for various values of V for diffusive 
coupling (single molecule). The lower right hand panel shows PLV as a function of V for all four types 
of coupling mechanisms. 


time period T by, 


P{t) 


1 

T 




t-T 


(14) 


The range of PLV value is [0,1]. When the value of PLV is zero, the two systems are uncoupled, whereas 
if the value of PLV is one then the two systems are perfectly phase synchronized [18]. 


Application: N. crassa circadian model 

The circadian model we consider is the simplified reduced model |43ll47] in which the negative feedback 
mechanism is incorporated during genetic regulation of the clock protein [26] as shown in Fig.l. We briefly 
describe the reaction mechanisms of the model as follows. The biochemical network of the circadian 
rhythm involve transcription and transportation of mRNA (M) with rate Vs by clock gene (frq) into 
cytosol. Then M decays with rate Vm or transported into cytosolic protein {Fc) with rate kg. Fc either 
decays with rate Vd or get inside the nucleus to form TV with rate fci which is a reversible reaction. This 
gives rise six reaction steps: (j) M, M (j), ^ Fq, Fc ^ (j), Fc ^ Fjv, Fn ^ Fc; with transition 

rates given by, xi = X 2 = X 3 = kgM, Xi = VdV Xs = kiFc and 

X6 = k 2 FM- The reactions are derived from the set of classical differential equations, which describe the 
molecular mechanisms, to stochastic reaction steps using transformation equation, r(V) = Fqdy)) where 
Pi and Xi are stochastic and classical transition rates, X = [M, Fc, and x = ^ = [^i fc, are 

state vector of the molecular populations and concentrations, and V is dimensionless system size. 

In the reduced reaction network model of N. crassa which exhibit circadian rhythm, the stochastic state 
of the system at any instant of time t can be described by a state vector X{t) = [M{t), Fc{t), Fjv(t)]“^. 
At any stochastic state of the system, the participating molecules in the network suffer decay or creation 
of molecules 120]. The time evolution of the transitional probability of the stochastic states is described 
by a Master equation I3135]. The stochastic dynamics was simplified by Gillespie based on identification 
of reaction number and reaction time, and we use this algorithm to simulate time evolution of molecular 
populations. 
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In our numerical simulation, we consider two dimensional array of oscillators {N x L) coupled by 
nearest neighbour diffusive or mean-field coupling as shown in Fig 2. We considered fixed boundary 
conditions i.e. for molecular species Fjv; = 0, = 0, = 0, = 0 and similarly 

the same condition is applied for the remaining other molecular species. So for any oscillator, there are 
at least two and at most four nearest neighbour coupled oscillators. 

Results and discussions 

The number of nuclei which can be viewed as self-sustained oscillators present in a single cytoplasm of 
N. crassa is small and finite. The oscillators are coupled via various coupling mechanisms mentioned in 
the previous section and the oscillators naturally prefer to choose the coupling mechanism that enable to 
process information quicker and easier i.e. in another word the mechanism which enable the oscillators 
synchronize fastest. In our large scale simulation we use stochastic simulation algorithm due to Gillespie 
m which we developed in Java language and we take 50 oscillators (nuclei) in a N. crassa single cell. 
During our simulation we also assume that these nuclei are static relative to each other because of the 
reasons that the rate of diffusion of the diffusing molecules or proteins (M, Fq etc) from one oscillator 
to another is much much faster as compare to the oscillator motion and there are other cellular processes 
such as microtubule, spindle etc that cause resistance in their relative motion. 

We first present the results of nearest neighbour single molecule bidirectional diffusive coupling among 
the 50 oscillators where M is taken to be coupling molecule as shown in Fig. 3. The simulation is 
done at system size V = 100 with coupling constant c = 1.4 and coupling is switch on at time = 
400hours. The upper two panels show the F^r dynamics as a function of time for 10 oscillators and 
their corresponding phase difference for each pairs of oscillators as a function of time calculated 
using Hilbert transform (|T^ . When the oscillators are uncoupled, the Fn dynamics of the oscillators 
and their corresponding Acj) are evolved independently of each other, whereas when the oscillators are 
coupled, the F^r dynamics exhibit the same correlated behaviour and their corresponding A(j) fluctuate 
about a constant value separating synchronized and desynchronized regimes. This phase transition like 
behaviour separating synchronized and desynchronized regimes is again verified by lowermost two panels 
the phase locking values, PLV dynamics of the oscillators and thier synchronization manifold recurrence 
plot. In desynchronized regime, the PLV evolved independently and in synchronized regime it remain at 
a constant value near the value 1. In the case of synchronized manifold recurrence plot, the synchronized 
regime the points are concentrated along the diagonal, whereas in the desynchronized regime, the points 
spread away from the diagonal. 

Now we compare the rate of synchrony for four different coupling mechanisms at constant system 
size, V = 100 and coupling constant c = 0.7 in Fig. 4. The results due to one molecular species diffusive 
coupling (coupling molecule is M) are shown in Fig. 2 such that (a) phase plot i.e. A(j> vs t, (b) phase 
locking value plot i.e. PLV vs t and (c) synchronized manifold i.e. recurrence plot. The plots show that 
the rate of synchrony is very weak because of the large fluctuation in the A(j) and PLV vs t plots in 
synchronized regime as well as larger spreading of points away from the diagonal. Similarly, the results 
for double molecular species diffusive coupling, single and double molecular species mean-field coupling 
mechanisms are shown in Fig. 4 (d), (e), (f); Fig. 4 (g), (h), (i); and Fig. 4 (j), (k), (1) respectively. 
The plots in all cases of coupling mechanisms show that the rate of synchrony is strongest in the case of 
double molecular species mean-field coupling and weakest in the case of single molecular species diffusive 
coupling. Our simulation results support our theoretical claim in the expression (fT^ . 

Next we present our simulation results for the four coupling mechanisms at constant system size 
V = 100 and for different values of coupling constants as shown in Fig. 3. The results of single molecular 
species diffusive coupling are shown in Fig. 5 (a), (b) and (c) respectively. It shows that for lower values 
of c, the A(/) and PLV evolve randomly as the function of time even if the coupling is switch on as shown 
in Fig. 5 (a) and (b). However this random fluctuation starts co-ordinating around a constant value as 
c increases and we found the rate of synchrony to be strongest at c = 1.4 and remain the same rate for 



Figure 7. Plot of phase diagram in which V for all four coupling mechanisms are shown as a function 
of their corresponding coupling constants c (c could be c, c', e or e' according to their respective 
coupling mechanisms) separating synchronized and desynchronized regimes. 


c > 1.4. This claim is supported by our recurrence plot in Fig. 5 (c). The results of double molecular 
species diffusive coupling, single and double molecular species mean-field coupling mechanisms are shown 
in Fig. 5 (d), (e), (f); Fig. 5 (g), (h), (i); and Fig. 5 (j), (k), (1) respectively. The values of coupling 
constants at which rate of synchrony is strongest in these three coupling mechanisms are found to be 
c = 1.1, 1.3 and 0.8 respectively. 

We show our simulation results for the single molecular species diffusive coupling mechanism at 
constant coupling constant c = 0.6 and for different values of system sizes in three panels of Fig. 6. 
The results show that even if the coupling is switch on at t = 400 hours, the dynamics of Ac/) and FLV 
randomly evolve with as a function of t for system sizes V (400. However the rate of synchrony is strongest 
at y = 400 and remain the same rate for system sizes V > 400. Our claim is again supported by the 
recurrence plot shown in third panel of Fig. 6. Then we compare the PLV values for different values of 
V for four different coupling mechanisms as shown in fourth panel of Fig. 6. The plots show that the 
rate of synchrony is strongest for double molecular species mean-held coupling scheme and weakest for 
single molecular species diffusive coupling scheme again supporting our theoretical claim (TT^ . 

We hnally compare the four coupling mechanisms in the phase diagram shown in Fig. 7 separating 
desynchronized and synchronized regime. The plot clearly indicates that the rate of synchrony in double 
molecular species mean-held coupling scheme is strongest and quickest since the area bounded by the 
curve in synchrony regime is largest. Whereas the rate of synchrony in single molecular species diffusive 
coupling scheme is weakest and slowest since the area bounded by the curve in synchrony regime is least. 
This clearly support our theoretical claim in (1121) again. 

Conclusion 

Naturally the coupled oscillators process information from one oscillator to another by selecting the 
coupling mechanism which gives fastest information processing and easily available in nature out of 
different coupling schemes. This idea of selection of coupling mechanisms by the coupled oscillators 
can be seen by measuring the rate of synchrony where this rate is maximum. In such situation of 
maximum rate of synchrony, the rate of information processing among the coupled oscillators is quickest 
and probably the agents which are responsible for information transfer among the oscillators may be 
easily available. 

In our large scale simulation for different coupling mechanisms that might happen in real situation of 
N. crassa circadian model, we found that mean-field coupling mechanism might be a prefered coupling 
scheme out of four schemes we studied. This numerical results support our theoretical claim also. The 
coupling schemes we studied so far allow relay or long range information transfer among the coupled 
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oscillators. However there are different coupling mechanisms such as direct coupling, delay time coupling 
etc which we do not study here thinking that either these coupling mechanisms are not realistic or slow 
coupling mechanisms specially in our N. crassa system. 
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